\documentclass{article}
\usepackage{lmodern}
\usepackage[a4paper,scale=0.8]{geometry}
\usepackage{listings}
\usepackage[dvipsnames]{xcolor}
\lstset{
    language=C++,
    basicstyle=\ttfamily,
    breaklines=true,
    keywordstyle=\bfseries\color{NavyBlue}, 
    commentstyle=\itshape\color{black!50!white},
    stringstyle=\bfseries\color{PineGreen!90!black},
    columns=flexible,
    numbers=left,
    numbersep=2em,
    numberstyle=\footnotesize,
    frame=single,
    framesep=1em
}

\usepackage{pgf-umlcd}

\usepackage{ctex}
\usepackage{amsmath}
\usepackage{graphicx}
\usepackage{float}
\usepackage{datetime}

\title{数值分析作业3}
\author{陈乐瑶  3210103710}
\date{\today}

\begin{document}
\maketitle

\part{样条的设计}

\begin{tikzpicture}
    \begin{class}[text width=9cm]{Spline}{0,0}
        \operation[0]{+solve(): void}
        \operation[0]{+operator()(const double \&point): double}
    \end{class}
    \begin{class}[text width=7cm]{ppformSpline}{-5,-3.5}
        \inherit{Spline}
        \attribute{+n: int}
        \attribute{+x: vector$<$double$>$}
        \attribute{+y: vector$<$double$>$}
        \attribute{+bct: BCType}
        \attribute{+polys: vector$<$Polynomial$>$}
        \attribute{+cond\_left: const double}
        \attribute{+cond\_right: const double}
        \operation{+ppformSpline(\_x,\_y,\_cond\_left,\_cond\_right)}
        \operation{+solve(): void}
        \operation{+operator()(const double \&point): double}
    \end{class}
    \begin{class}[text width=6cm]{Bspline}{5,-3.5}
        \inherit{Spline}
        \attribute{+n: int}
        \attribute{+x: vector$<$double$>$}
        \attribute{+y: vector$<$double$>$}
        \attribute{+x\_complete: vector$<$double$>$}
        \attribute{+bct: BCType}
        \attribute{+bases: vector$<$BsplineBase$>$}
        \attribute{+cond\_left: const double}
        \attribute{+cond\_right: const double}
        \operation{+Bspline(\_x,\_y,\_cond\_left,\_cond\_right)}
        \operation{+solve(): void}
        \operation{+operator()(const double \&point): double}
    \end{class}
    \begin{class}[text width=8cm]{Function}{0,-17}
        \operation[0]{+operator()(const double \&x): double}
        \operation{+get\_diff(const double \&x): double}
        \operation{+get\_secdiff(const double \&x): double}
    \end{class}
    \begin{class}[text width=7cm]{Polynomial}{-5,-12}
        \inherit{Function}
        \attribute{+n: int}
        \attribute{+coef: vector$<$double$>$}
        \operation{+Polynomial()}
        \operation{+operator()(const double \&x): double}
        \operation{$\ldots$}
    \end{class}
    \begin{class}[text width=7cm]{BsplineBase}{5,-12}
        \inherit{Function}
        \attribute{+N: int}
        \attribute{+knots: vector$<$double$>$}
        \operation{+BsplineBase(\_N,\_knots)}
        \operation{+operator()(const double \&x): double}
        \operation{$\ldots$}
    \end{class}
    \unidirectionalAssociation{ppformSpline}{}{}{Polynomial}
    \unidirectionalAssociation{Bspline}{}{}{BsplineBase}
\end{tikzpicture}

\section{Function类}
\begin{itemize}
    \item 纯虚函数operator ()，用于计算函数值，每个子类都需要实现。
    \item 虚函数get\_diff ()和get\_secdiff ()，用于计算函数的一阶和二阶导数，默认由中心差分实现，子类可以选择重载。
\end{itemize}

\subsection{Polynomial类 —— Function的子类}
\begin{itemize}
    \item 用于表示多项式函数$p(x) = c_0 + c_1x + c_2x^2 + \cdots + c_nx^n$，实现分段多项式样条插值。
    \item 包含成员变量n和coef，分别表示多项式的阶数和系数。
    \item 实现基本运算的重载。
\end{itemize}

\subsection{BsplineBase类 —— Function的子类}
\begin{itemize}
    \item 用于表示B样条基函数$B_{i,k}(x)$。
    \item 包含成员变量N和knots，分别表示基函数的阶数和节点向量。
    \item 实现operator ()时，对$N = 1$的情况直接处理，对$N > 1$的情况递归调用。
\end{itemize}

\section{Spline类}
\subsection{ppSpline类 —— Spline的子类}
\begin{itemize}
    \item 用于实现分段多项式样条插值方法。
    \item 用模板类实现一次和三次样条函数，其中三次样条函数实现Complete，Second和Natural三种边界条件，Natural边界条件
    只需设置左右边界条件默认值为0。
    \item 用多项式类的数组储存每一段的多项式函数。
    \item 样条插值时使用了第一章的Hermite插值方法，对于三次样条函数，需要计算每一段的一阶导数，因此需要重载Function类的
    get\_diff ()函数。
\end{itemize}

\subsection{三次ppSpline类的设计原理——以Complete边界条件为例}
由（3.3）（3.4）和（3.15）式得到线性方程组：
\begin{equation}
    \begin{bmatrix}
        2 & \mu_2 & 0 & \cdots & 0 \\
        \lambda_3 & 2 & \mu_3 & \cdots & 0 \\
        0 & \lambda_4 & 2 & \cdots & 0 \\
        \vdots & \vdots & \vdots & \ddots & \vdots \\
        0 & 0 & 0 & \cdots & 2
    \end{bmatrix}
    \begin{bmatrix}
        m_2 \\
        m_3 \\
        m_4 \\
        \vdots \\
        m_{N-1}
    \end{bmatrix}
    =
    \begin{bmatrix}
        3\mu_2 f[x_2, x_3] + 3\lambda_2 f[x_1, x_2] - \lambda_2 m_1 \\
        3\mu_3 f[x_3, x_4] + 3\lambda_3 f[x_2, x_3] \\
        3\mu_4 f[x_4, x_5] + 3\lambda_4 f[x_3, x_4] \\
        \vdots \\
        3\mu_{N-1} f[x_{N-1}, x_N] + 3\lambda_{N-1} f[x_{N-2}, x_{N-1}] - \mu_{N-1} m_N
    \end{bmatrix}
\end{equation}
其中$m_1$和$m_N$为边界条件，$\lambda_i = \frac{x_i - x_{i-1}}{x_{i+1} - x_{i-1}}, \mu_i = \frac{x_{i+1} - x_i}{x_{i+1} - x_{i-1}}, i = 2, 3, \ldots, N-1$。\\
下面以三次Complete边界条件pp样条函数为例，给出部分代码实现：
\begin{lstlisting}
    if (bct == Complete) {
        SparseMatrix<double> A(n - 2, n - 2);
        VectorXd b(n - 2), r(n - 2);
        vector<Triplet<double>> tripletlist;

        //fill the matrix A and middle part of vector b
        for (int i = 0; i < n - 2; i++) {
            tripletlist.push_back(Triplet<double>(i, i, 2));
            if (i != 0) {
                double lambda = (x[i + 2] - x[i + 1]) / (x[i + 2] - x[i]);
                tripletlist.push_back(Triplet<double>(i, i - 1, lambda));
            }
            if (i != n-3) {
                double mu = (x[i + 1] - x[i]) / (x[i + 2] - x[i]);
                tripletlist.push_back(Triplet<double>(i, i + 1, mu));
            }
            b(i) = 3 * (x[i + 1] - x[i]) / (x[i + 2] - x[i]) * (y[i + 2] - y[i + 1]) / (x[i + 2] - x[i + 1]) + 3 * (x[i + 2] - x[i + 1]) / (x[i + 2] - x[i]) * (y[i + 1] - y[i]) / (x[i + 1] - x[i]);
        }
        A.setFromTriplets(tripletlist.begin(), tripletlist.end());
        A.makeCompressed();

        //fill the first and last row of vector b
        b(0) -= (x[2] - x[1]) / (x[2] - x[0]) * cond_left;
        b(n - 3) -= (x[n - 2] - x[n - 3]) / (x[n - 1] - x[n - 3]) * cond_right;
        
        SparseLU<SparseMatrix<double>> solver;
        solver.compute(A);
        r = solver.solve(b);

        //calculate the coefficients with the result of linear equation
        VectorXd m(n);
        m(0) = cond_left;
        m(n - 1) = cond_right;
        m.segment(1, n - 2) = r;
        for (int i = 0; i < n - 1; i++) {
            vector<double> x_h = { x[i], x[i], x[i + 1], x[i + 1] };
            vector<vector<double>> y_h = { { y[i], y[i], y[i + 1], y[i + 1] }, { m[i], m[i], m[i + 1], m[i + 1] } };
            Hermite herm(x_h, y_h);
            polys.push_back(herm.get_poly());
        }
    }
\end{lstlisting}


\subsection{Bspline类 —— Spline的子类}
\begin{itemize}
    \item 用于实现B样条插值方法（包括Cardinal B样条插值）。
    \item 用模板类实现一次、二次和三次B样条函数，其中二次B样条函数默认为Theorem3.58中的Cardinal B-Spline设计，
    三次B样条函数则实现Complete，Second和Natural三种边界条件，Natural边界条件只需设置左右边界条件默认值为0。
    \item 用B样条基函数类的数组储存支持区间基函数。
\end{itemize}

\subsection{三次Bspline类的设计原理——以Complete边界条件为例}
由Theorem3.49和Lemma3.27，$f(x_i) = a_{i-2}B_{i-2}^3(x_i) + a_{i-1}B_{i-1}^3(x_i) + a_i B_{i}^3(x_i) $。\\
又由边界条件和B样条基函数的导数性质，$f'(x_1) = a_{-1}(\frac{3B_{-1}^2(x_1)}{x_1-x_{-2}} - \frac{3B_{0}^2(x_1)}{x_2-x_{-1}}) + a_0(\frac{3B_{0}^2(x_1)}{x_2-x_{-1}} - \frac{3B_{1}^2(x_1)}{x_3-x_{0}}) + a_1(\frac{3B_{1}^2(x_1)}{x_3-x_{0}} - \frac{3B_{2}^2(x_1)}{x_4-x_{1}})$，
$B_{-1}^2(x_1) = B_{2}^2(x_1) = 0$，$f'(x_n) = a_{n-3}(\frac{3B_{n-3}^2(x_n)}{x_{n-2}-x_{n-5}} - \frac{3B_{n-2}^2(x_n)}{x_{n-1}-x_{n-4}}) + a_{n-2}(\frac{3B_{n-2}^2(x_n)}{x_{n-1}-x_{n-4}} - \frac{3B_{n-1}^2(x_n)}{x_{n}-x_{n-3}}) + a_{n-1}(\frac{3B_{n-1}^2(x_n)}{x_{n}-x_{n-3}} - \frac{3B_{n}^2(x_n)}{x_{n+1}-x_{n-2}})$，$B_{n-1}^2(x_n) = B_{n}^2(x_n) = 0$。\\
由此可得线性方程组：
\begin{equation}
    \begin{bmatrix}
         - \frac{3B_{0}^2(x_1)}{x_2-x_{-1}} & \frac{3B_{0}^2(x_1)}{x_2-x_{-1}} - \frac{3B_{1}^2(x_1)}{x_3-x_{0}} & \frac{3B_{1}^2(x_1)}{x_3-x_{0}} & \cdots & 0 \\
        B_{-1}^3(x_1) & B_0^3(x_1) & B_1^3(x_1) & \cdots & 0 \\
        0 & B_0^3(x_2) & B_1^3(x_2) & \cdots & 0 \\
        \vdots & \vdots & \vdots & \ddots & \vdots \\
        0 & 0 & 0 & \cdots & - \frac{3B_{n}^2(x_n)}{x_{n+2}-x_{n-1}}
    \end{bmatrix}
    \begin{bmatrix}
        a_{-1} \\
        a_0 \\
        a_1 \\
        \vdots \\
        a_n
    \end{bmatrix}
    =
    \begin{bmatrix}
        f'(x_1) \\
        f(x_1) \\
        f(x_2) \\
        \vdots \\
        f'(x_{n})
    \end{bmatrix}
\end{equation}

下面以三次Complete边界条件Bspline样条函数为例，给出部分代码实现：
\begin{lstlisting}
    if (bct == Complete) {
        SparseMatrix<double> A(n + 2, n + 2);
        VectorXd b(n + 2), r(n + 2);
        vector<Triplet<double>> tripletlist;

        //fill the middle rows
        for (int i = 1; i < n + 1; i++) {
            for (int j = i - 1; j < i + 2; j++) {
                double tmp = bases[j](x_complete[i + 2]);
                tripletlist.push_back(Triplet<double>(i, j, tmp));
            }
        }

        //fill the first and last row
        vector<BsplineBase> bases2;
        for (int i = 0; i < n + 1; i++) {
            vector<double> knots2(x_complete.begin() + i+1, x_complete.begin() + i+5);
            BsplineBase base2(2, knots2);
            bases2.push_back(base2);
        }
                
        double tmp1 = -3 * bases2[0](x_complete[3]) / (x_complete[4] - x_complete[1]);
        double tmp2 = 3 * bases2[1](x_complete[3]) / (x_complete[5] - x_complete[2]);
        double tmp3 = -3 * bases2[n - 1](x_complete[n + 2]) / (x_complete[n + 3] - x_complete[n]);
        double tmp4 = 3 * bases2[n](x_complete[n + 2]) / (x_complete[n + 4] - x_complete[n + 1]);
        tripletlist.push_back(Triplet<double>(0, 0, tmp1));
        tripletlist.push_back(Triplet<double>(0, 1, - tmp1 - tmp2));
        tripletlist.push_back(Triplet<double>(0, 2, tmp2));
        tripletlist.push_back(Triplet<double>(n + 1, n - 1, tmp3));
        tripletlist.push_back(Triplet<double>(n + 1, n, - tmp3 - tmp4));
        tripletlist.push_back(Triplet<double>(n + 1, n + 1, tmp4));

        A.setFromTriplets(tripletlist.begin(), tripletlist.end());
        A.makeCompressed();

        //fill the vector b
        b[0] = cond_left;
        b[n + 1] = cond_right;
        for (int i = 1; i < n + 1; i++) {
            b[i] = y[i - 1];
        }

        SparseLU<SparseMatrix<double>> solver;
        solver.compute(A);
        coef = solver.solve(b);
    }
\end{lstlisting}


\subsection{CurveFitting类 —— 样条插值的高维补充}
\begin{itemize}
    \item 用于实现用一次和三次的两种样条函数进行高维曲线拟合的方法。
    \item 用length表示Definition 3.71中的cumulative chordal length，用于计算每一段的参数值。
    \item 用knots表示高维节点，用curve表示高维拟合点。
    \item 用curve\_fit ()函数表示为每一维拟合一条以cumulative chordal length为独立参数的样条曲线。
    实现输入一组高维点，对每一维根据要求进行样条插值，得到多组拟合点，并将点储存在curve中。
    \item 其中三次样条的边界条件默认为Natural，便于在离散函数不易求导的情况下进行拟合。
\end{itemize}

\subsection{三次CurveFitting类的设计原理}
由Definition3.71和Algorithm3.72，对于每一维，计算每一段的参数值，然后根据参数值计算每一段的系数，最后得到每一维的样条函数。\\
下面以三次pp样条函数为例，给出部分代码实现（其中边界条件默认为Natural）：
\begin{lstlisting}
    for (int i = 0; i < dim; i++) {

        //x is the knots of each dimension
        vector<double> x;
        for (int j = 0; j < n; j++) {
            x.push_back(knots[j][i]);
        }
        if (splinetype == ppform) {
            if (N == 3) {
                ppformSpline<3, Natural> ppspline_3(length, x);

                //nodes is the fitted points of each dimension
                vector<double> nodes;
                for (double k = 0; k <= length[n - 1]; k += dl) {
                    nodes.push_back(ppspline_3(k));
                }

                //curve is the set of fitted points
                curve.push_back(nodes);
            }
        }
    }
\end{lstlisting}


\part{样条的测试}
\section{A：$f(x) = \frac{1}{1+25x^2}$}
用两种方法对$f(x) = \frac{1}{1+25x^2}$进行三次样条插值，在不同插值点个数N的条件下，插值结果如下：
\begin{figure}[htbp]
    \begin{minipage}[t]{0.3\textwidth}
        \includegraphics[width=5.5cm]{./pic/A1_6.png}
        \caption{N = 6}
    \end{minipage}
    \begin{minipage}[t]{0.3\textwidth}
        \includegraphics[width=5.5cm]{./pic/A1_11.png}
        \caption{N = 11}
    \end{minipage}
    \begin{minipage}[t]{0.3\textwidth}
        \includegraphics[width=5.5cm]{./pic/A1_21.png}
        \caption{N = 21}
    \end{minipage}
    \begin{minipage}[t]{0.35\textwidth}
        \includegraphics[width=6cm]{./pic/A1_41.png}
        \caption{N = 41}
    \end{minipage}
    \begin{minipage}[t]{0.35\textwidth}
        \includegraphics[width=6cm]{./pic/A1_81.png}
        \caption{N = 81}
    \end{minipage}
    \caption{cubic spline interpolation of $f(x) = \frac{1}{1+25x^2}$}
\end{figure}

max-norm error随插值点个数的变化如下：
\begin{figure}[htbp]
    \begin{minipage}[t]{0.4\textwidth}
        \includegraphics[width=6cm]{./pic/A_error.png}
        \caption{误差曲线}
    \end{minipage}
    \begin{minipage}[t]{0.4\textwidth}
        \includegraphics[width=6cm]{./pic/A_error_output.png}
        \caption{误差输出}
    \end{minipage}
\end{figure}

由Theorem3.12的误差分析，$|f(x)-s(x)| \leq \frac{1}{16}h^4\max_{x\in[a,b]}|f^{(4)}(x)|$，其中$h = \max_{i=1}^{N-1}(x_{i+1}-x_i)$。
因此我们令$h = \frac{1}{N-1}$，作$\log(h)$和$\log(error)$的图像如下：
\begin{figure}[htbp]
    \includegraphics[width=6cm]{./pic/A_error_convergence.png}
    \caption{log (error) of cubic spline interpolation of $f(x) = \frac{1}{1+25x^2}$}
\end{figure}
由图像可以看出斜率为4，符合误差分析的结果，验证了样条插值的误差收敛率。

\section{C：$f(x) = \frac{1}{1+x^2}$}

用二次和三次Cardinal B样条对$f(x) = \frac{1}{1+x^2}$进行插值，插值图像如下：
\begin{figure}[htbp]
    \includegraphics[width=6cm]{./pic/C.png}
    \caption{Cardinal B-spline interpolation of $f(x) = \frac{1}{1+x^2}$}
\end{figure}

在D中给定的点集上，误差曲线和误差输出如下：
\begin{figure}[htbp]
    \begin{minipage}[t]{0.4\textwidth}
        \includegraphics[width=6cm]{./pic/C_error.png}
        \caption{误差曲线}
    \end{minipage}
    \begin{minipage}[t]{0.4\textwidth}
        \includegraphics[width=6cm]{./pic/C_error_output.png}
        \caption{误差输出}
    \end{minipage}
\end{figure}

问题：\\
Q1：为什么部分点的误差接近机器精度？\\
A1：因为-3.5，-0.5，0.5，3.5四点为二次基数B样条的插值点，而-3，0，3三点为三次基数B样条的插值点，因此这些点的误差接近机器精度。\\
Q2：哪一种样条插值方法的更精确？\\
A2：三次样条插值方法更精确，由误差曲线可以明显看出，且由误差分析也可知，因为三次样条插值方法的误差收敛率为4，而二次样条插值方法的误差收敛率为3。

\section{E：heart function}
用三次pp样条对heart function进行插值，因为heart function在[0,1]上的导数不易求，因此用Natural边界条件进行插值，在不同插值点个数n的条件下，插值结果如下：\\
\begin{figure}[htbp]
    \begin{minipage}[t]{0.33\textwidth}
        \includegraphics[width=5.5cm]{./pic/E_10.png}
        \caption{n = 10}
    \end{minipage}
    \begin{minipage}[t]{0.33\textwidth}
        \includegraphics[width=5.5cm]{./pic/E_40.png}
        \caption{n = 40}
    \end{minipage}
    \begin{minipage}[t]{0.33\textwidth}
        \includegraphics[width=5.5cm]{./pic/E_160.png}
        \caption{n = 160}
    \end{minipage}
    \caption{curve fitting of heart function}
\end{figure}
由图像可以看出，在给定范围内随n的增大，拟合效果越来越好。

\section{鸣谢}
本次作业的完成像在考试周对自己的追杀，特此鸣谢以下人员救我于水火：\\
感谢助教对本次作业提供的环境配置以及对C++问题的耐心解答，让我不断了解自己智商的下限。\\
感谢室友对本次作业的讨论和帮助，让我在段错误中没有自我了结并帮我删了好几版代码。\\
特别感谢coco老师，也就是copilot老师，让我在一次次tab键中感受gpt的魅力。

\end{document}